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Abstract 

Oh ' We present explicit filtration/backprojection-type formulae for the inversion of the 

spherical (circular) mean transform with the centers lying on the boundary of some 
polyhedra (or polygons, in 2D). The formulae are derived using the double layer po- 
tentials for the wave equation, for the domains with certain symmetries. The formulae 
are valid for a rectangle and certain triangles in 2D, and for a cuboid, certain right 
prisms and a certain pyramid in 3D. All the present inversion formulae yield exact 
reconstruction within the domain surrounded by the acquisition surface even in the 
' presence of exterior sources. 

oo ■ 

OO 

Introduction 

ON 

The spherical mean Radon transform and related inversion formulae are important, in partic- 
ular, for solving problems of thermo- and photo-acoustic tomography (TAT/PAT) [23, 24, 28]. 
In these modalities a region of interest (ROI) within a human body is subjected to a very 
short electromagnetic (EM) pulse which causes thermoelastic expansion of the tissues and 
generates an outgoing acoustic wave. The resulting acoustic pressure is measured by the 
detectors placed on a surface surrounding ROI. The initial pressure distribution f(x) is not 
uniform: tumors absorb much more EM energy than healthy tissue and thus generate a 
much stronger signal. Hence, by recovering f(x) one obtains valuable medical information. 

Under certain simplifying assumptions the problem of finding the initial pressure can 
be re-stated as that of inverting the spherical mean Radon transform in 3D over spheres 
with centers lying at the measuring surface. A similar problem in 2D arises when instead of 
usual point-like detectors one measures the acoustic pressure with linear detectors [6, 31, 32]. 
A detailed discussion of the solvability, range conditions, and inversion techniques for the 
spherical mean transform in 2D and 3D can be found in [1-4,11-13,18-21,33,35-37] and 
references therein. 

In the present paper we derive explicit inversion formulae for the spherical (or, in 2D, cir- 
cular) mean transform for certain measurement surfaces with special symmetries. In general, 
such explicit formulae play an important role in inverse problems: they provide an important 
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theoretical insight and, in some cases, they serve as a starting point for the development of 
efficient reconstruction algorithms. For example, a very popular filtration/backprojection 
(FBP) algorithm (see, for example, [25]) is obtained as a result of discretization of one of 
the explicit inversion formulae for the 2D Radon transform. 

For the spherical mean Radon transform the inversion formulae of the filtration/back- 
projection type are known only for special types of the acquisition surfaces. In particular, 
several different formulae for the spherical measuring surface have been found in [38, 40] (in 
3D), in [11] (in odd dimensions), in [10] (in even dimensions), and in [21] (in arbitrary dimen- 
sions). A general family containing as particular cases all the above-mentioned formulae for 
a sphere, was derived in [26]. The so-called "universal" backprojection formula [40] (already 
mentioned above) holds not only for a spherical measuring surface, but also for a the surface 
of an infinite cylinder and for a plane [39,40] (both in 3D). Other formulae for the planar 
acquisition are also known [6,8,9,27,29]. The existence of the FBP-type formulae for these 
three types of surfaces (and the absence of such formulae for other smooth surfaces) can be 
explained by certain symmetries common to a sphere, a plane, and a cylinder, as discussed 
in [20]. 

For other measuring surfaces there exist certain numerical techniques for solving the 
problems TAT/PAT (see [1, 20, 33, 35-37] and references therein); but the inversion formulae 
for such surfaces are not known. (The notable exception is the explicit general solution 
obtained in [1] for any closed surface; however, this solution is given in terms of functions 
of the wave operator rather than in a form of a closed form integro-differential expression). 
In the present paper we derive explicit FBP-type inversion formulae for the spherical mean 
transform with the centers lying on the surfaces of certain polyhedra (in 3D) or on the 
boundaries of certain polygons (for circular mean transform in 2D). 

The paper is organized as follows. In Section 1 we provide a formal description of the 
problem and supply some introductory information. Next, in Section 2 we explicitly solve 
the time reversal problem for the wave equation in a cube or a square (in 2D) by means of 
the double layer wave potentials. Such a solution is possible due to certain symmetries of 
these regions. Next, in Section 3, from the explicit solution of the wave equation we derive 
FBP-type formulae for the inversion of the spherical (circular) mean transform with centers 
on the boundary of a cube (a square). The present inversion formulae have an interesting 
property of being insensitive to acoustic sources lying outside of the region surrounded by the 
acquisition surface, as discussed in Section 3.3. Numerical examples of image reconstructions 
obtained using our formulae are provided in Section 3.4. Finally, in Section 4 we apply the 
same methods to some other domains possessing the necessary symmetries. In 2D these 
include a rectangle and three types of triangles: an equilateral and a right isosceles triangles, 
and a triangle with the angles |, | and ^. In 3D our formulae extend to cuboids, right 
prisms whose base is one of the three triangles mentioned above, and a pyramid whose side 
faces are equal right isosceles triangles. 

1 Formulation of the problem 

Throughout the paper we assume that the source of the acoustic pressure is supported within 
a bounded region fl, the detectors are placed on the boundary dfl, and the speed of sound is 
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constant and equal to 1 in the whole space M. n . The time-dependent acoustic pressure p(x,t) 
satisfies the wave equation 

Ptt = A x p, t>0, xeW 1 
p(x,0) = f(x), Pt (x,0) = 0, 

where initial pressure f(x) is supported within Q. The detectors measure the values P(y,t) 
of the pressure on the boundary: 

P(y,t)=p(y,t)\ ye9a , te[0,oo). (2) 

Our goal is to reconstruct f(x) from measurements P(y,t). 

When the measurements are done using conventional point-like acoustic detectors, the 
above problem has to be solved in 3D setting. However, when the so-called line detectors 
are utilized [6,31,32], the 3D problem reduces to a set of similar 2D inverse problems for 
the wave equation, so from the practical point of view both 2D and 3D cases are of interest. 

In the 3D case the measurements need only be done during the time interval [0, dianrsl], 
since the signal vanishes after t = diamf2 due to the Huygens principle [34]. In the 2D 
case, the values of P(y,t) in the whole interval t G [0, oo) are needed for the theoretically 
exact reconstruction, and we will assume that these values are known. (However, if P(y, t) 
is known for all t G [0, dianrsl], the remaining values for the time interval [diamQ, oo) can be 
explicitly reconstructed, see [30] and the discussion at the end of the next section). 



1.1 Spherical means 

The problem of finding f(x) can also be re- formulated in terms of inverting the spherical 
mean transform. Indeed, in 3D the solution to the initial value problem (1) can be written 
in the form of the Kirchhoff formula 

p(x,t) = -^tM(x,t), (3) 

where M(x,r) are the spherical means of f(x) 

M(x,r) — — J f(x + ru)du), 
s 2 

and where S 2 is the unit sphere in 3D. Thus, using (2) one can relate the spherical means 
M(y,t) to P(y,t): 

P(y,t) =-tM(y,t), 
and, taking into account that M(y, 0) = 0, solve for M(y, t) 

M{y,t) = \ f P(y,t')dt'. 
1 Jo 

Now the problem is reduced to finding function f(x) from its spherical means M(y,t) with 
centers on dQ. 
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m(x,r) — — I f(x + r(jj)du. 



Similarly, in 2D, one can introduce the circular means m(x,r) : 

1 

2^ 

Now solution of (1) can be expressed through m(x,r) by the formula 

t 

m(x, r)rdr 



s 1 



Vt 2 - r 2 



and, in particular, 



o 

The integral in the right hand side of the latter equation is the well known Abel transform 
whose inversion formula is also well known. Thus, circular means m(y,r) can be found from 
P(y,t) by the following formula 



r 

m(y,r) = ^- J 



o 



P(y,t)dt 



which again leads to the problem of reconstructing f(x) from its cylindrical means m(y,r). 
(Since m(y,r) vanish for r > diamf2, the knowledge of P(y,t) in the time interval t G 
[0, diamfi] is sufficient to reconstruct m(y,r) using the latter formula; moreover, values of 
P(y,t) for t > diamf2 can now be recovered from m(y,r) by equation (4).) 

1.2 Reconstruction strategies 

It is well known that the initial pressure f(x) can be found by the time reversal, i.e. by 
solving the wave equation back in time. In particular, in the 3D case, since solution to (1) 
vanishes within Q after t = diamf2 due to the Huygens principle, one can solve the following 
initial-boundary- value problem backwards in time : 

'u tt = A x u, te[0,T\, xettcR 3 , 
u(x,T) = 0, ut(x,T) = 0, 
u(y,t) = P(y,t), yedtl, 
T = diamf2. 

Function f(x) we seek to recover equals u at the moment 0: 

f(x)=u(y,0). (6) 

In 2D the solution to the wave equation does not vanish in a finite time, and therefore 
the initial condition should be replaced by the condition of the decrease at t — > oo : 

u tt = A^m, t e [0, oo), xeftc K 2 , 

lim u(x,T) = 0, lim u t (x,T) = 0, (7) 

T^oo T-5>oo v ' 

u(y,t) = P(y,t), ye da 
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As in 3D, f(x) is found using formula (6). 

The above initial-boundary- value problems can be solved numerically using, for example, 
finite difference methods (e.g. [5,7]). For certain acquisition surfaces methods based on 
eigenfunction expansion can also be very effective [22]. Our goal, however, is to obtain 
an explicit filtration/backprojection type formula for the solution. If the existing detectors 
could measure (in addition to the pressure) the normal derivative of the pressure, such an 
explicit analytic solution could have been represented, using the Green's formula, in the 
form of the double- and single- layer potentials. Unfortunately, the normal derivative is not 
known. Still, explicit representations of the solution of the wave equation by means of the 
double layer potential are possible for surfaces with certain symmetries, such as, for example, 
a sphere, a cylinder, and a plane. In fact, the so-called "universal backprojection formula" 
[40] for these surfaces can be interpreted as evaluation of a certain double layer potential. 

In the next section we find explicit inversion formulae for the case when Q is a cube 
(in 3D) or a square (in 2D). First, we obtain explicit solutions for the problems (5) and 
(7) in the form of the double layer potentials. These formulae allow us to reconstruct f(x) 
from the measured values of the pressure P(y,t); in Section 3 they also will be used to 
derive explicit filtration/backprojection time formulae for recovering f(x) from its spherical 
(circular) means in these domains. 



2 Explicit solution of the wave equation in a cube or a 
square 

2.1 Double layer potentials 

In 3D, the retarded free space Green's function G 3D+ (x, t) is given by the following for- 
mula [34] 

Air\x\ 

it solves the equation 

G 3 t D+ - A X G 3D+ = 6(t)6(\x\) (8) 

subject to the radiation condition as t — > oo. We will need the advanced Green's function 
G 3D ~(x,t), 

G™-(x,t) = G 3D ^x,-t) = S -^±^, (9) 



that solves (8) back in time and satisfies the radiation condition as t — > — oo. 

The double- and single- layer potentials are frequently used in the theory of partial 
differential equations (see e.g. [34]). The wave potentials that solve the time dependent 
wave equation, were introduced in [14, 15] (see also [16]). We summarize below some of the 
basic properties of the wave potentials that will be needed for further use. (Here we sacrifice 
generality for simplicity and discuss only the properties that will be used later in the text). 

Given a plane II and a function ip(y,t), y E U which is continuous, twice differentiable in 
time and compactly supported in both variables, we introduce the (advanced) single layer 
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potential U(x,t) and the double-layer potential V(x,t) with the density (p(y,t) as follows: 



U(x,t) 



J J \f(y,t')G 3D -(y-x,t-t')ds(y)dt', (10) 



i 1 n 



V(x,t) = J J ip ( y ,t')^G 3D -(y~x,t-t')ds(y)dt', (11) 
r 1 n v 

where n is the normal to II, and ds(y) is the standard area element. Note that both U(x,t) 
and V(x,t) satisfy the wave equation for all x outside II; both potentials describe waves 
propagating (backwards in time) away from the support of ip, and they satisfy the radiation 
condition at infinity as t — >■ — oo. Also, V(x,t) = for all x e II, and due to the standard 
jump conditions [16] 

Inn V(y±en,t) = ±]-(p(y,t), yell. 

Since G 3D ~(x, t) vanishes for positive values of t, equation (11) implies that if tp(y, t) vanishes 
for t > t then so does V(x, t). Finally, since n is constant, using (9) V(x, t) can be re-written 
as follows 

v(x , t) = _i . v /fO^+kzi!)^). (12) 

An J 4:71 \y — x\ 

n 

= -j-div fn ^^y-^ dsiy). 
4tx J 4Tj\y-x\ yy> 
n 

We notice that, due to (12), (for a planar surface only), 

V(x,t) = -n-V x U(x,t). (13) 
In 2D, the advanced free space Green's function G 2D ~(x,t) has the following form [34]: 

r, \x\ < \t\, t < 0, 



G 2D -(x,t) = | 2,V^-' 



otherwise 



Further, for a given straight line L and for a continuous, twice different iable in time, com- 
pactly supported in both variables function <p(y, t), y G L we define the single-layer potential 
U(x,t) with the density (p(y,t) : 

oo 

U(x,t) = J J \p(y,t')G 2D (y-x,t-t')dt'dl(y) 

L t 

oo 

hi l ^' t + T) ^- (x -y/ Tdl{l/) 

^^-(x-yY'dTdHy). (14) 



L \ x -y\ 



± f f fdtp(y,r) 
2tt J J \dr r 

L \x-y\ 



r=t+r 
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Now, instead of dealing directly with -§^G (x,t) (which would lead to strongly singular 
integrals), we define the double layer potential V(x,t) by means of (13) and (14): 

V(x,t) = -n ■ V x U{x,t) 



= —n ■ V x 
2tt 



L \ x -y\ 



9 <p(y,r) 
dr r 



a/t 2 - (x - y) 2 drdl(y) 



r=t+r 



1 

"2tt 



L \x-y\ 



1 

"2tt 



L \x-y\ 



9 <p(y,r) 

dr r 

9 <p(y,r) 
dr r 



n 



V y y/T*-(x-y)*dTdl(y) 



r=t+T 



n-(x-y) 



-.drdl(y), 



where n is the normal to L, and dl(y) is the standard arc length. As before, V(x,t) satisfies 
the wave equation (now in 2D) for all x outside L and it describes an outgoing wave prop- 
agating (backwards in time) away from the support of ip. Also, V(x,t) = for all x G L, 
and 

Inn V(y±en,t)=±-(p(y,t), y G L. 

e— >0,e>0 Z 



2.2 Time reversal in a cube 

An explicit solution to the initial-boundary-value problem (5) for the case when Q is a cube 
can be obtained using the double layer potentials, as described below. 

Without loss of generality let us assume that Q = (—a, a) x (—a, a) x (—a, a), and 



boundary <9f2 consists of six faces Sj, j 



,6: 



on = |J Sj. 

In this case the wave leaves f2 after time T = diamfi = 2\/3a. Let us assume that x = 
(xi,X2,xs), and that the face Si lies in the plane x\ = a. First, we will solve the initial- 
boundary-value problem with non-zero boundary conditions on face Si only: 



'u tt = A x u, t G [0, T], x G Q, C M 3 , 
u(x,T) = 0, u t (x,T) = 0, 
lim £ _> ,£>o u(y - en u t) = P(y,t), y G Si, 

6 

u(y,t) = 0, y G U Sj, 

j=2 



(15) 



where rii is the exterior normal to Si. The solution procedure we are about to outline will 
apply as well to the other faces Sj, j = 2, 6, and the solution of problem (5) will then be 
readily obtained as the sum of the solutions corresponding to each face. 
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To simplify the presentation let us introduce an operator E that extends any function 
h(x 2 ,x 3 ) defined in the open square (—a, a) x (—a, a) to a function h Tcpl (x 2 , x 3 ) defined in 
R 2 by means of odd reflections. In detail, if h repl = E 2D h, then h rcpl coincides with h within 
(—a, a) x (—a, a); further, for any point (2:2, £3) lying in an open square s n>m with the side 
2a and centered at (2na,2ma), h Tcp[ (x 2 ,x 3 ) is defined by the following formula: 

h Tep \x 2 ,x 3 ) = {-l) m+n h{a, {-l) n {x 2 - 2nd), {-l) m {x 3 - 2ma)). 

Finally, on the boundary of each square s n ^ m we define function h mpl (x 2 ,x 3 ) to be zero. 

Let us denote by II" the plane containing Si, and introduce plane IT^ parallel to II" 
and passing through the point (—3a, 0,0). We equip both planes with the same normal 
rti coinciding with the exterior normal to face Si. We introduce function P™ pl (x 2 , x 3 , t) = 
E 2D (P(x,t)\ Si ) , i.e. we replicate the boundary values of u on face Si. Now, for t G [0,T], 

we define the double-layer potentials V®{x,t) and V*(x,t) with the density — 2P 1 rcpl (y 2 , IJ3, t) 
supported on planes II" and IlJ: 

V?(x,t) = ^div f |(a>y ^p^ Pi epl (y 2 ,y3^ + I (a, y 2 , y 3 ) - x\)dy 2 dy 3 , 
Vf{x,t) = £div/ 

R 2 

Note that since P[ ep[ is finitely supported in [0, T] in the time variable, the integrands in the 
above expressions vanish for large values of y 2 or y 3 , and therefore the integration is actually 
performed over bounded subsets of M 2 lying within the ball B(x,T — t). In order to make 
expressions (16) simpler, we will abuse notation by writing P[ cxA {y,t) for y G II" or y G Tl\ 
instead of P[ cv \y 2 ,y 3 ,t) : 

V 1 a ' b (x,t) = ^div I -^-Pr\y,t+\y-x\)ds{y), (17) 

2-k J \y-x\ 

Ul' b nB(x,T-t) 

where ds(y) is the standard area element, and V^' h means either V® or V\ b . 

The combined potential Vi(x,t) = V£{x,t) + V±{x,t) has some interesting properties. 
Let us denote by S 2 the face of cube VL opposite to Si. We notice that planes II" and IlJ 
are symmetric with respect to S 2 , and, due to the choice of the normals, potential Vi(x,t) 
vanishes on S 2 . Next, due to the oddness of P[ epl , potential Vi(x,t) vanishes on the other 
four faces Sj of the cube (j = 3, 6). Finally, since T is smaller than the distance between 
II" and Tl\, potential V^(x,t) equals on II" for all t <T. Therefore, 

lim V 1 {y-en 1 ,t)= lim V?{y - en x , t) = P(y, t), y G S ± . 

£— »0,£>0 £— »0,£>0 

Thus, potential Vi(x,t) solves initial-boundary- value problem (15). For brevity, we will 
re- write the expression for Vi (x, t) in the following form 

Vi(x,t) = ^div / j^—Pl^iy, t+\y- x\)ds(y). 

2n J \y — x\ 

(n°un5)ns(x,T-t) 
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In order to obtain the desired solution to (5), we replicate the boundary values on the 
faces S2, Sq : 

P/ epl (-,-,t)=^(p(-,t)| Sj ),j = 2,...,6. (18) 

Next, we introduce planes II", 11^ with normals rij, and define double layer potentials Vj with 
the densities -2P- cp \y, t) on IT?, U): 

Vj (x,t) = ±dW [ ^L^P^\ y , t +\y-x\)ds(y), j=2,...,6. 

Z7T J \y — x\ 

(njun5)ns(x,T-t) 

The sum of these potentials V(x, t) = Y^=i Vj( x > solves the initial-boundary- value problem 
(5), and, therefore, f(x) equals to V(x,0) : 

/(x) = X)^0) = ^ F divX; / n /^^p^ daCy). (19) 
i=1 J=1 (n?un5)ns(x,r) 

Interestingly, while each of the potentials V J a,b (x, t) solves the wave equation in the sense 

6 

of distributions (due to the discontinuities in the edges of Q) the sum ^2 Vj(x,0) of these 

3=1 

potentials solves the wave equation in the classical sense (if f(x) G Cq(Q)). Thus, we have 
proven the following 

Theorem 1 The Cq(Q) initial condition f(x) of the initial-boundary-value problem (1) in 
the case n = 3 can be reconstructed from the boundary values P(y,t) (defined by (2)) by 
formulae (19) and (18). 

Formula (19) closely resembles the so-called "universal backprojection formula" [40] es- 
pecially the version applicable to the measurements done over an infinite plane. However, in 
the present case the measurements are performed over a bounded surface (that of the cube), 
and the exact reconstruction at any point x is obtained by integration over a bounded set 

£(z,r)nU- =1 (n?uiTj). 



2.3 Time reversal in a square 

A similar technique can be used to find the explicit solution to the initial-boundary-value 
problem (7) in the case when Q is a square. The geometry of the problem is actually easier, 
thanks to the lower dimensionality of the problem. However, in 2D the Huygens principle 
does not hold and the solution will result from the integration over unbounded sets. 

We will assume that Q is a square (—a, a) x (—a, a), Si is the side of Q corresponding to 
Xi = a, and S2 is the side contained by the line Xi = —a. 

First, let us find solution u(x,t) to the following initial-boundary- value problem: 
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u tt = A x u, t G [0, T], xGOc K 2 , 
w(x,T)=0, u t (x,T) = 0, 
lim u(y-en,t) = P(y,t)ri(T,t), y E S u (20) 

[u(s/,t) = 0, yeS,-, J = 2,3,4, 

where r)(T,t) is a C°° cut-off function equal 1 for t G [0, T — 1] and vanishing with all the 
derivatives at t — T. We introduce operator E lD that extends a function h(x) defined in the 
interval (—a, a), to a function h repl (x) defined on R 1 by means of odd reflections. In detail, 
within the interval (—a, a) function /i repl coincides with h. For any x ^ (2m+l)a, m G Z, there 
is a unique integer such that |x — 2na\ < a. We define h rcpl (x) = (—l) n ^ h((—l) n ( x \x — 
2n(x)a)). Finally, at the odd integer points h rep[ vanishes: /i repl ((2m + l)a) = 0, Vm G Z. 

Let us call Li the straight line containing S\. Introduce the family of lines $i consisting 
of all straight lines parallel to Li, and passing through the points (a + 4ka, 0), fc G Z. 
(Obviously, Li G $1). The exterior normal rii to the side Si of the square will be used as 
the normal to all the lines in $1. Let us define function P 1 repl (x 2 ,t) by replicating values of 
P(x, t) for x G Si : 

P^W) = E 1D (P((a,x 2 ),t)). 

We will abuse notation by writing Pl epl (x, t) instead of the more accurate P 1 repl (x2(a;), t). Fur- 
ther, let us introduce the double layer potential Vi(x, t, T) with the density — 2P 1 repl (a;2, t)r)(T, t) 
supported on all lines in the family $1, defined by the formula 



00 , 

V 1 (x,t,T) = lj J ( 



dr r 



*i \x-v\ 



^■^-y) dTdm . 



r=t+T V t2 -( x ~ y) 2 



Note that in the above formula the integration in y is actually performed over a bounded 
subset of $1 due to the finite support of r\ in the second variable. 

Potential Vi(x,t,T) satisfies the wave equation in Q. Since p* cpl is odd with respect to 
the straight lines containing all sides of dQ, V^(x,t,T) vanishes on dfl. However, Vi(x,t,T) 
is not continuous across Si. The jump is such that 

lim Vi(y-em,t,T) = P^{y,t)ri(T,t) = P(y,t) V (T,t), y G S u t G [0.T]. 

£— »0,£>0 

Therefore, Vi(x,t,T) solves the problem (20). 

The next step is to obtain a solution u T (x,t) of the following boundary problem: 

'uf t = A x u T , te[0,T], xencM 2 , 

u T (x,T)=0, uf(x,T) = 0, (21) 
u T (y,t)=P(y,t) V (T,t), yedQ. 

This can be achieved by defining functions Pj Cpl (-,t) = E 1D ^P(-,t)\ s ^j and families of 
parallel lines $j with normals rij, j = 2,3,4, and by introducing the corresponding double 
layer potentials Vj(x, t, T), in a fashion similar to our definition of P[ epl , $1, rii, and Vi(x, t, T). 
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Now each Vj(x,t,T) satisfies the boundary problem with the inhomogeneous conditions on 
the side j, and the homogeneous conditions on the other three sides of the square. In turn, 
the sum of these potentials u T 



( 



4 4 

u T (x, t) = J2 *> T ) = " div E 

3 I j=l 



T-t 



\ 



n,— ; al{y)dT 



\ 



o 



V 7 " 2 - (y - x ) 2 



\y—x\ <t 



J 



7T ^ 



d Ff i {y,r)r } {T,r ) 
dr r 



*j |x-j/| 



r=i+r 



■ (a; - |j 
a/t 2 - (x - y) 2 



drdl{y) 



solves (21). 

It is known [17] that, as T grows to infinity, u T (x, 0) converges uniformly to /(x): 



/(x) = lim -u T (x 



4 00 / 

.«)=^E//( 



<9 ^ repl (y, r)?7(r, r) \ rij -(x-y) 



dr 



A/r 2 -{x-y) 



9 Pp P \y,r)\ n r (x-y) 



dr 



®j \x-y\ 



^/r 2 -(x-y) 



-drdl(y) 



zdrdl(y) 
(22) 
(23) 



The latter formula is very similar to the explicit inversion formula for measurements done 
from an infinite line [6]. Formula (23) can also be re-written using the Green's function 
G 2D ~(x — y,r) to illustrate its connection to the surface potentials: 

fW = lfl J jn r (x-y) (^ PjP( r y,r) ) G 2D ~(x-y,r)drdl(y). 

Theorem 2 The Cq(Q) initial condition f(x) of the initial-boundary-value problem (1) in 
the case n = 2 can be reconstructed from the boundary values P(y,t) (defined by (2)) by 
formulae (23). 



3 Explicit inversion of the spherical means: a cube and 
a square 

3.1 Inversion formula in 3D 

In 3D, in order to reconstruct f(x) from its spherical means M(y,t) with centers on dQ, we 
introduce replicated spherical means M* epl (y,t) defined as follows 

Mf pl (-,-,t) = E 2D (M(.,t)\ Sj ) ,j = 1,...,6. (24) 
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Since the replication operator E 2D is invariant with respect to the time variable, from (3) 
we conclude that 

Pr P \y>t)=§- t tM^\y,t). (25) 

Now formula for reconstructing f(x) from spherical averages in 3D is obtained by substituting 
(25) into (19): 

J 1 (n?un5)ns(x,T) 
We thus have proven 

Theorem 3 In 3D, a Cq(Q) function f(x) supported within the cube Q can be reconstructed 
from its spherical means M(y,t) with centers on dQ by the formula (26), where Mj Cpl (y,t) 
is defined by (24)- 

Formula (26) has a form of the filtration by differentiation, followed by the backprojection, 
followed by the divergence operator. It is similar to the inversion formula obtained by the 
author [21] for the inversion of the spherical means with centers on a sphere (in 3D). In the 
latter case the integration is over the surface of a sphere, and no replication of the data is 
needed. 



/ n 



ld_ 

tdt 



tM^\y,t) 



ds(y). 



(26) 



t=\cr.— n\ 



3.2 Inversion formula in 2D 

In this section we derive an explicit inversion formula for finding f(x) from the values of 
its circular means with the centers lying on the boundary of a square. The starting point 
for the derivation is equation (23) that reconstructs the sought function from the boundary 
values of the solution of the wave equation; we re-write it in the following form 



4 °c 



i =1 $ 

4 



dPr\y,r)\ 1 



dr 



\Jr 2 — s 2 



dr 



dl(y) 



s=\x—y\ 



(27) 



where 



d P; cp \y,r ) 

dr r 



\Jr 2 — s 2 



dr. 



(28) 



First, we will express Pj epl (y,s) in terms of the circular averages m(y,r). Let us apply the 
extension operator to m(y,r) and define rn J epl (y,r) as follows: 



m?*(-,r) = E 1D (m(-,r)\ Sj ), J = l,...,4. 



(29) 
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Function P rcpl (-,t) is related to m repl (-,r) in the same way as P(-,t) is related to m(-,r): 



d r m™ p \y,r)rdr d 
~dt J y/W^ 2 ' 



dt 



\Jt 2 — r 2 — m^ cpl (y, r)dr 



t^mf p \y,r) 



dr 



Vt 2 - r 2 



rdr 



t 



where we integrated by parts twice, and used the fact that ra^ pl {y, 0) = for all y. Now one 
can compute the derivative needed in (28): 



t 



d f d ( d 



-m cp \y, r))Vt 2 - r 2 dr = t \ 
at J or \ror J J 

o o 



t d_ ( _d 

dr \ rd 



-mJ p \y,r)J 



dr. 



By substituting the above equation into (28), interchanging the integration order, and notic- 
ing that mj epl (y,r) = for r > 2y/2a, one obtains the following expression for P™ pl (y, s) : 



t t 



d ( d 



s 
2%/2e 



t 



= lim 



did rcpl 

./ Wr\rTr m > M 
o 

2\p2a 



T 



drdt 



t 



max(s,r) 



Vt 2 - rWt 2 - s 2 



di 



dr 



= lim 



ez. I Tr (Jk<^- r >) [ ln + VW ~^) ' \ ln (I 



r 2 - s 2 \) 



dr 



1 

2 



o 

2v / 2a 



ln(|r 2 -s 2 |) 



d ( d 



m Tepl (y, r) ) dr 



dr \ rdr 



+ lim / — 



d ( d 



t^oo / dr V rdr 3 



-m 



rcpl 



(y, r)^j ln (Vt 2 - s 2 + VT 2 - r 2 ) dr 



(30) 



The last integral in (30) vanishes: 

2v / 2a 



lim f 4- ( 4r m T P \y, r)\ In ( VT 2 - s 2 + VT 2 - r 2 ) dr 

t^oo J dr \ror 3 J \ J 



lim 



d 



m™ pl (y,r) 



t^oo J rdr 3 vy ' (VT 2 - s 2 + VT 2 - r 2 ) VT 2 - r 



-.dr = 0, 
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where we integrated by parts and used the finite support of m^ cp (y, r) in r. Now (30) can be 
re- written in the following form 

2\/2a 2\/2a 

P;j\y,s) = - 1 - J In (|r 2 - S 2 |) |-(^ m ^( y , r ))d r = PV / ^^ { y,r )d r, 

o o 

where "PV" indicates that the integral is understood in the principal value sense. Substitu- 
tion of the latter expression for Pp^iy, s) into (27) yields the following sequence of equivalent 
reconstruction formulae for f(x): 



/(*) 



j =1 i 



n,; 



[x-y) 



2^2a 



PV 



1 d 



dl(y) 



2ir ^ J \x 



(x - y) 



( 2\/2a 

I- / ln(|r 2 - s 2 | 



- l -t 



n.; 



ds 



\ 



d 



s=\x-y\ 
\ 



—mf pl (y,r)dr 



(x - y) 



\x 



, 

ds 



o 

2\/2a 



rm 



rcpl 



(y,r) 



dr 



J 



J 



dl(y) 



dl(y) 



s=\x—y\ 



s=\x—y\ 



2\p2a 



— div 

7T 



4 

E 



PV 



2\/2a 

PV 



rm 



rcpl 



(y,r) 



r 2 — s 2 



dr 



dl(y) 



m 



rcpl 



(y,r) 



o 



r 2 — (x — y) 



-rdr 



s=\x—y\ 



dl(y). 



(31) 



(32) 



(33) 



Thus, the following statement holds: 

Theorem 4 In 2D, a C 2 (f2) function f(x) supported within the square Q can be recon- 
structed from its circular means m(y,r) with centers on dfl by the formula (33) (or, equiv- 
alently, by formulae (31) or (32)), where m^ cpl (?/,r) is defined by (29). 

Formula (33) closely resembles formula 8.19 [2] that reconstructs a function from its 
circular means centered on a circle. 

The drawback of formulae (31)-(33) is that the integration has to be done over an un- 
bounded set Uj =1 Qj. However, by analyzing the expression in parentheses in (31) one can 
notice that the integrand of the outer integral decreases as \x — y\~ 3 for large values of y. 
Moreover, for values of \x — y\ > 2\[2a the integrand is infinitely smooth. Therefore, by in- 
tegrating over a subset $ trunc containing all points y G U^ =1 $j such that dist(y, Q) < diamf2 
one correctly reconstructs the singularities of f(x), and obtains a good quantitative approx- 
imation to f(x). This conclusion is supported by the numerical evidence presented in the 
next section. 
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3.3 Reconstruction in the presence of exterior sources 

The time reversal reconstruction methods in TAT/PAT (i.e. methods based on the numerical 
solution of the wave equation backwards in time) have the following interesting property. If 
the support of the initial perturbation f(x) has a part lying outside the region enclosed by 
the acquisition surface dfl, the time reversal still correctly reconstructs the part of f(x) lying 
within fi. In other words, the exterior sources do not affect the reconstruction within f2 [1, 22]. 
Since reconstruction methods based on the expansions in the series of the eigenfunctions of 
the Dirichlet Laplacian on Q [1, 22] are theoretically equivalent to the time reversal, they also 
inherit the insensitivity to the exterior sources. On the other hand, all previously known 
filtration/backprojection-type inversion formulae for the spherical mean transform do not 
have this property: in the presence of an exterior source the reconstruction within fl will be 
incorrect. 

The present inversion formulae (and their generalizations in the following sections) are 
based on the explicit solution of the wave equation (by means of double layer potentials). 
Thus, they inherit from the time reversal methods the insensitivity to the exterior sources. 
A numerical example demonstrating accurate reconstruction in the presence of an exterior 
source is presented in the next section. 

3.4 Numerical examples 

An example illustrating reconstruction of a function from its spherical means with centers 
lying on the surface of a cube is presented in Figure 1. As a phantom we used a set of 25 
characteristic functions of balls of various radii supported within the cube (—1, 1) x (—1, 1) x 
(—1, 1). The centers of all the balls lied in the plane x% = 0. The cross section of our phantom 
by the latter plane is shown in Figure 1(a). The part (b) of this figure demonstrates the cross 
section of the image reconstructed on the grid of size 129 x 129 x 129 using formula (26). The 
simulated detectors (i.e. the centers of the integration spheres) were placed on the uniform 
Cartesian 129 x 129 grids on each of the faces of the cube; there were 257 integration spheres 
for each detector with the radii varying uniformly from to 2\/3. Our algorithm was based 
on a straightforward discretization of equation (26), where finite differences were used to 
compute the derivatives. 

In order to illustrate the reconstruction in the presence of exterior sources we added to 
the phantom shown in Figure 1(a) a ball of radius 0.08 located at the position (1.1,0,0), 
and computed the spherical means corresponding to this extended set of sources. The size 
and location of this additional ball were chosen so that all the spherical means used for the 
reconstruction still vanished when the radii exceeded 2\/3. The result of the reconstruction 
computed using formula (26) is shown in Figure 1(c). The resulting image may seem indis- 
tinguishable from the one in Figure 1(b); in fact, the difference is about 4% in the L°° norm. 
(In the absence of the discretization error the images would exactly coincide). 

In Figure 2 we demonstrate the reconstruction of a function in 2D from its circular means 
centered on the boundary of a square (—1, 1) x (—1, 1) using a modified formula (31), with 
the integration restricted to the subset $ trunc of U| =1 $j contained within the disk of radius 
3\/2 centered at the origin. (All the points y satisfying the condition dist(y,fi) < diamfi 
are contained in $ trunc ). The main goal of this experiment was to evaluate the effect of 
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(a) (b) (c) 

Figure 1: Reconstruction in 3D from spherical means centered on a surface of a cube (a) 
phantom (b) reconstructed function (c) reconstruction in a presence of a source outside of 
the cube 




(a) (b) 

Figure 2: Reconstruction in 2D from circular means centered on the perimeter of a square 
(a) phantom (b) reconstruction 

truncating the integration region. To this end we took care to eliminate other sources 
of error. In particular, as a phantom we chose a smooth function shown in Figure 2(a). 
Moreover, our implementation of formula (31) used a spectrally accurate algorithm for 
computing the Hilbert transform and a higher order polynomial interpolation during the 
backprojection. (The Hilbert transform arises since fraction r/(r 2 —s 2 ) in 31 can be re- written 
as 0.5/(r — s) + 0.5/ (r + s).) The reconstructed function is shown in In Figure 2(b). The 
image is almost identical to the phantom; the relative reconstruction error in this example 
equals 7.4 • 10~ 3 in the L°° norm. In other words, the truncation error is negligible for most 
practical applications, if all the points y satisfying the condition dist(y, Q) < diamf2 are 
included in the truncated domain. 



16 



4 Explicit inversion formulae for some other domains 



4.1 Re-visiting the replication procedure for a cube and a square 

The inversion formulae introduced in this paper for a cube (in 3D) and a square (in 2D) 
are based on the explicit representations of a solution of the wave equation by double- 
layer potentials. In general, double layer potentials are frequently used to solve Dirichlet 
problems for various PDE's in the interior and/or exterior of bounded regions. Usually such 
formulations lead to boundary integral equations for the potential densities that have to be 
solved numerically. However, as explained in the previous sections, if the detectors (centers 
of spherical means) are located on the boundary of a cube (or, in 2D, a square), the solution 
of the wave equation can be represented by the double layer potentials supported on certain 
periodic planes, with certain periodic densities. Due to the symmetries in these geometries, 
the potentials vanish on the planes, and the jump conditions relate the densities directly to 
the known boundary values, thus providing the desired explicit reconstruction formulae. 

This approach can be extended to obtain explicit reconstruction formulae for the mea- 
surements made from the boundaries of certain polygons (in 2D) and polyhedra in (3D) 
possessing sufficient symmetries. These formulae will be similar to the present formulae 
for the square and the cube, but the planes supporting the potentials and the replication 
procedure defining the densities will be different in each case. 

Before extending the inversion formulae to these new domains, let us provide an alterna- 
tive description of the replication procedure for a square and a cube - such that it would be 
easier to generalize to other geometries. We start with a square, in 2D. Our goal is, again, 
to find the solution u T (x,t) of the initial-bouindary-value problem (21) provided the bound- 
ary values of u are known. To this end, let us construct a tesselation of the plane 1R 2 by 
squares of the size 2a x 2a, one of which coincide with Q = (—a, a) x (—a, a). In each square 
we define -u repl to equal ±u T (in some local coordinates). The signs are chosen so as to be 
different in any two adjacent squares. The local coordinates are defined in such a way that 
the function in any two adjacent squares is odd with respect to the joint side of these two 
squares. The resulting function -u repl (x, t) is continuous in x everywhere except the boundary 
of each square. At each point of the boundary it has a jump equal to a doubled negative 
limiting value of u repl (x, t) as x approaches the boundary from the inside. Now we define the 
double layer potential supported on all lines containing the boundary of the squares, with 
the density equal to the jump in -u rcpl . This new procedure leads to the same double layer 
potentials as were described in Section 2.3, with the densities equal —2Pj Cpl (x2,t)i](T,t) on 
each family of lines <3>j. The rest of the reconstructing procedure is the same. 

Similarly, in order to obtain the double layer potential for a solution of the wave equation 
in a cube Q = (—a, a) x (—a, a) x (—a, a), we tessellate the plane with shifted versions of Q. 
In each shifted cube we define u rcpl (x, t) as a reflected version of u(x, t) with a sign opposite 
to the sign used in the neighboring cubes. The local coordinates in each cube are defined in 
such a way that u repl (x, t) is odd with respect to each side of each cube. Now we define the 
double layer potentials supported on all planes that contain all faces of all cubes, with the 
density equal to the jump of u Tcpl (x, i) across the cube faces. This new procedure produces 
the same set of the double layer potentials as was utilized in the Section 2.2, and it leads to 
the same reconstruction formulae. Although the double layer potential described above is 
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(a) (b) (c) 

Figure 3: The tesselation schemes for some triangular domains (a) equilateral triangle (b) 
right isosceles triangle (c) triangle with angles |, | and | 

supported on an unbounded set, only a bounded subset of it is actually used to reconstruct 
the image within fl, due to the finite speed of sound. 

One obvious extension of the present inversion formulae is to rectangular domains in 2D 
and to cuboids in 3D. All the symmetries needed for the double layer potentials to form 
an explicit solution remain in place in these cases, and all the reconstruction formulae (for 
the wave equation and for the spherical/circular means) remain the same, with a proper 
re-definition of integration lines and planes. 

4.2 Inversion formulae for certain triangular domains in 2D 

In 2D, the inversion formulae presented in Sections 2.3 and 3.2 can be generalized to certain 
triangular domains. Namely, a equilateral triangle, a right isosceles triangle, and a triangle 
with the angles |, | and | can be used to tessellate the plane while preserving the necessary 
symmetries. The corresponding tesselations are shown in Figure 3. In each of the three 
cases shown in the figure, the sides of the triangular fundamental domain Q are marked by 
different types of lines (solid, dashed, and dotted ones) so that the orientation of the local 
coordinates in each of the triangles is clear from the picture. Signs "+" and "-" indicate 
the choice of u or — u in each triangle. Clearly, the resulting u repl (x,t) is odd with respect 
to all lines containing all sides of all the triangles. Therefore, if one defines the double layer 
potential with the density supported on these lines and equal to the jump of u repl (x, t) across 
the corresponding face, the value of this potential will vanish on all the triangle sides, and, 
due to the jump condition, this potential will solve the initial-boundary- value problem for the 
wave equation in Q. Therefore, the analog of the formula (23) is valid with the corresponding 
re-definition of p re P' and with Q being one of the above mentioned triangles. 

This, in turn, allows us to obtain the analog of the inversion formulae (31) - (33) for 
the reconstruction from circular means, by defining m repl in a similar fashion. As it is the 
case for the square domain, the latter formulae are exact if the integration is done over the 
unbounded graph. However, due to the fast decrease of the expression in the parentheses 
in (31) one can hope that the integration over the region containing all points y satisfying 
the condition dist (2/, SI) < diamfi will produce an accurate approximation to f(x). 
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Figure 4: Extending u(x,t) by odd reflections to a cube 

4.3 Inversion formulae for certain polyhedral domains in 3D 

In addition to cuboids, the techniques of the previous sections can also be used to find explicit 
solutions and inversion formulae in 3D for some other domains. In particular, one can derive 
explicit inversion formulae for reconstructing a function from its spherical means centered 
on the surface of a right prism of height H, whose base is one of the triangles discussed in 
the previous section. As before, it may be easier to start with finding the explicit solution 
to the wave equation u(x, t) whose values on the boundary of the prism are known. In order 
to obtain such a solution, we replicate u(x,t) in the directions parallel to the base of the 
prism, according to the schemes in Figure 3. Then we replicate the resulting slab of height 
H in the odd fashion with respect to the planes containing the bases of the prism, until 
u(x,t) is extended to function u repl (x,t) defined in the whole M. 3 in x. Further, we define 
the double layer potential with the density supported on planes containing the faces of all 
the replicated prisms, and equal to the jumps of u(x, t) across the replicated faces. This 
double layer potential vanishes on all the faces, and, due to the jump condition, provides a 
solution to the initial-boundary- value problem for the wave equation in 3D. This results in 
an explicit formula analogous to (19). In order to obtain the analog of the formula (26) it is 
enough to define M repl according to the replication procedure described above. Due to the 
Huygens principle, the integration is actually performed over the subset of the replicated 
faces satisfying the condition \x — y\ < diamfi; both formulae (for the spherical means and 
for the wave equation) are exact in this case. 

Finally, the present method can be extended to the triangular pyramid Q whose vertices 
have coordinates (0, 0, 0), (a, 0, 0), (0, a, 0), and (0,0, a). (This region can also be described 
as a pyramid whose side faces are equal right isosceles triangles.) Again, starting with the 
solution u(x,t) of the wave equation in Q, we define u repl (x, t) as follows. First, we reflect 
u(x,t) in odd fashion with respect to the vertical faces of Q to extend u to a rectangular 
pyramid (see Figure 4). Then we reflect (oddly) the pyramid with respect to its triangular 
faces to obtain six pyramids (only four are shown in the figure) to form a cube. Then the cube 
is used to tessellate the whole space. As before, we define the double layer potential by the 
jumps of the resulting function u rcpl (x,t) across all the replicated faces. The symmetries of 
such an object guarantee that the potential vanishes on all the faces, and the jump conditions 
yield the required limiting values as one approaches the boundary of Q from the inside. Thus, 
one obtains the exact solution to the initial-boundary-value problem corresponding to the 
time reversal of the wave equation. By implementing a similar replication procedure to define 
M repl , one obtains an analog of the inversion formula (26) for reconstructing a function from 
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its spherical means centered on the surface of f2. 

It is worth noting that, similarly to the formulae for a cube and a square, all the formulae 
discussed in the present and previous sections also have the property of being insensitive to 
sources lying outside of the domain surrounded by the acquisition surface. 

Conclusions 

By utilizing the double layer wave potentials we found explicit solutions to the initial- 
boundary-value problems for the wave equation in certain 2D and 3D domains. In prob- 
lems of TAT/PAT this yields explicit expressions for the result of the time reversal, and 
thus explicitly reconstructs the initial pressure distribution f(x). Further, by formulating 
the problem in terms of spherical (circular) means, we obtain from the solutions to the wave 
equations explicit closed- form inversion formulae of filtration/backprojection type for the 
corresponding domains. In 2D, we presented inversion formulae for a square (or a rectan- 
gle), and we outlined their generalization for an equilateral triangle, right isosceles triangle, 
and a triangle with the angles |, | and |. In 3D, we derived formulae for a cube (with im- 
mediate extension to a cuboid), and we discussed their generalization to right prisms whose 
bases are one of the three above-mentioned triangles, and to a pyramid whose side faces are 
equal right isosceles triangles. 

In all cases the formulae for the inversion of the wave equation data are in the form of 
the double layer potentials supported on a series of planes, with the densities whose values 
equal to the known boundary values replicated according to the procedures described in the 
text. The formulae for inverting the spherical (or circular) means easily follow from the ones 
for the wave data. In the 3D case, due to the Huygens principle the integration in all the 
formulae is restricted to a bounded subset of the planes. In 2D, in order to obtain the exact 
reconstruction one has to integrate over an unbounded domain. However, our numerical 
experiments suggest that for a sufficiently accurate reconstruction it is enough to integrate 
over a bounded subset containing all points y satisfying the condition dist(y, Q) < diamfi. 

Finally, all the presented inversion formulae remain exact within the domain even in the 
presence of sources located outside the domain. (In order to observe this phenomenon in 
3D one may need to extend the measurement time (or the radii of the spherical means) to 
make sure that the signal vanishes after the observation time). Interestingly, although the 
time reversal and some series methods also have this property, the present formulae are the 
only known explicit FBP-type inversion formulae that exhibit such behavior. 
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